!--------------------------------
! Programa dedicado a calcular
! derivadas ateh terceiro grau
! atraves de expansoes de taylor.
! A precisao eh relatada em arquivo
! que sai.

! MODELO DE ARQUIVO DE ENTRADA AO FINAL
! DO CODIGO

! As seguintes derivadas são calculadas:
! 
! - Derivada para frente de 2 pontos
! df = (f_1 - f_0) / h
!
! - Derivada para tras de 2 pontos:
! df = (f_0 - f_-1) / h
!
! - Derivada simetrica de 3 pontos:
! df = (f_1 - f_-1) / 2h
!
! - Derivada simetrica de 5 pontos:
! df = (f_2 -8*f_1 +8*f_-1 - f_-2) / (12*h)
!
! -Derivada segunda simetrica de 3 pontos:
! ddf = (f_1 - 2*f_0 + f_-1) / (h^2)
! 
! - Derivada segunda simétrica de 5 pontos
! ddf = (-f_2 + 16*f_1 -30*f_0 +16*f_-1 -f_-2) / (12*h^2)
!
! - Derivada terceira anti-simetrica de 5 pontos
! dddf = (f_2 - f_1 + f_-1 - f+-2) / (2*h**3)

! A seguinte funcao foi utilizada:
! f = exp(4x)cos(x/2)
! No ponto x = 1/3


! Programa compilado e testado com:
! $ gfortran -o exerA exerA.f90
! $ ./exerA

! Forma de utilizacao auto-explicativa
! atravez de comentarios.

! _o_o_ oOo _o_o_ 23/10/2010


PROGRAM exerA
  IMPLICIT NONE
  INTEGER :: N,  i, ids3, idf2, idt2, ids5, idds5, idddas5
  REAL (KIND=8) :: ds3, df2, dt2, ds5, dds5, dddas5, df, ddf, dddf, x=1./3
  REAL (KIND=8) :: bds3=1E10, bdf2=1E10, bdt2=1E10, bds5=1E10, bdds5=1E10, bdddas5=1E10
  REAL (KIND=8), DIMENSION(20) :: H
  REAL (KIND=8), DIMENSION(-2:2) :: F

  ! Calculando as derivadas no ponto x
  df=exp(4*x)*( 4*cos(x/2) - sin(x/2)/2)
  ddf=exp(4*x)*((63./4)*cos(x/2)-4*sin(x/2))
  dddf=exp(4*x)*(61*cos(x/2)-(191./8.)*sin(x/2))

  ! Lendo os dados
  ! 1) Numero de intervalos
  OPEN(11,file='tab1_in.dat')
  read(11,*) N

  ! 2) Intervalos em si
  i=1
  DO WHILE (i<=N)
    read(11,*) H(i)
    i=i+1
  END DO

  open(unit=10,file="tab1_out.dat")
  write(10,*) "h ", "|sim 3pt ", "|frente 2pt ", "|tras 2pt ", "|sim 5pt ", "|seg sim 5pt ", "|ter anti-sim 5pt"

  ! Calculando as derivadas atraves das
  ! aproximacoes e dos erros.
  i=1
  F(0)=exp(4*x)*cos(x/2)
  write(*,*) "iniciando loop"
  DO WHILE (i<=N)
    ! Calculando o valor das funcoes nos pontos em questao
    F(-2)=exp( 4*(x-2*H(i)) )*cos( (x-2*H(i))/2)
    F(-1)=exp( 4*(x-H(i)) )*cos( (x-H(i))/2)
    F(1)=exp( 4*(x+H(i)) )*cos( (x+H(i))/2)
    F(2)=exp( 4*(x+2*H(i)) )*cos( (x+2*H(i))/2)


    ! Calculando as derivadas, erros e registrando as menores ocorrencias.
    ! Veja no cabecalho para detalhamento das formulas

    ! Derivada simetrica de 3 pontos
    ds3=(F(1)-F(-1))/(2*H(i)) - df
    IF (ds3 < bds3) THEN
      bds3=ds3
      ids3=i
    END IF

    ! Derivada para frente de 2 pontos
    df2=(F(1)-F(0))/H(i) - df
    IF (bdf2 > df2) THEN
      bdf2 = df2
      idf2=i
    END IF

    ! Derivada para tras de 2 pontos
    dt2=(F(0)-F(-1))/H(i) - df
    IF (abs(bdt2) > abs(dt2)) THEN
      bdt2 = dt2
      idt2=i
    END IF


    ! Derivada simetrica de 5 pontos
    ds5=-(F(2)-8*F(1)+8*F(-1)-F(-2))/(12*H(i)) -df
    IF (abs(bds5) > abs(ds5)) THEN
      bds5 = ds5
      ids5=i
    END IF

    ! Derivada segunda simetrica de 5 pontos
    dds5=(-F(2)+16*F(1)-30*F(0)+16*F(-1)-F(-2))/(12*H(i)**2) -ddf
    IF (abs(bdds5) > abs(dds5)) THEN
      bdds5 = dds5
      idds5=i
    END IF

    ! Derivada terceira anti-simetrica de 5 pontos
    dddas5=(F(2)-2*F(1)+2*F(-1)-F(-2))/(2*(H(i)**3)) - dddf
    IF (abs(bdddas5) > abs(dddas5)) THEN
      bdddas5 = dddas5
      idddas5=i
    END IF

    i=i+1
  END DO

  ! Mostrando no terminal qual a melhor escolha e porque
  write(*,*) "---------------------------------"
  write(*,*) "---------------------------------"
  write(*,*) "Para derivada simetrica de 3 pontos,"
  write(*,*) "o melhor valor de h eh", H(ids3)
  write(*,*) "Pois o erro eh: ", bds3

  write(*,*) "Para derivada para frente em 2 pontos,"
  write(*,*) "o melhor valor de h eh", H(idf2)
  write(*,*) "Pois o erro eh: ", bdf2

  write(*,*) "Para derivada para tras em 2 pontos,"
  write(*,*) "o melhor valor de h eh", H(idt2)
  write(*,*) "Pois o erro eh: ", bdt2

  write(*,*) "Para derivada simetrica em 5 pontos,"
  write(*,*) "o melhor valor de h eh", H(ids5)
  write(*,*) "Pois o erro eh: ", bds5

  write(*,*) "Para derivada segunda simetrica em 5 pontos,"
  write(*,*) "o melhor valor de h eh", H(idds5)
  write(*,*) "Pois o erro eh: ", bdds5

  write(*,*) "Para derivada terceira anti-simetrica em 5 pontos,"
  write(*,*) "o melhor valor de h eh", H(idddas5)
  write(*,*) "Pois o erro eh: ", bdddas5
end PROGRAM

!!!!!!!!!!!!!!!!!!
! MODELO DE ARQUIVO DE ENTRADA.
! DESCOMENTE e COPIE em tab1_in.dat
! 14
! 0.5
! 0.2
! 0.1
! 0.05
! 0.01
! 0.005
! 0.001
! 0.0005
! 0.0001
! 0.00005
! 0.00001
! 0.000001
! 0.0000001
! 0.00000001
